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Abstract 

We  derive  a  general  expression  that  quantifies  the  total  entanglement  production  rate  in  continuous 
variable  systems,  where  a  source  emits  two  entangled  Gaussian  beams  with  arbitrary  correlators.  This 
expression  is  especially  useful  for  situations  where  the  source  emits  an  arbitrary  frequency  spectrum, 
e.g.  when  cavities  are  involved.  To  exemplify  its  meaning  and  potential,  we  apply  it  to  a  four-mode 
optomechanical  setup  that  enables  the  simultaneous  up-  and  down-conversion  of  photons  from  a 
drive  laser  into  entangled  photon  pairs.  This  setup  is  efficient  in  that  both  the  drive  and  the 
optomechanical  up-  and  down- conversion  can  be  fully  resonant. 


1.  Introduction 


Entanglement  is  an  essential  feature  of  quantum  mechanics  and  a  crucial  resource  for  quantum  communication 
and  information  processing.  The  most  common  situation  involves  a  source  that  continuously  produces 
entangled  beams.  One  of  the  most  natural  characteristics  of  such  a  source  is  obviously  the  rate  at  which  it 
generates  entanglement.  If  the  source  sends  out  pairs  of  entangled  particles,  with  subsequent  pairs  completely 
independent,  this  rate  can  simply  be  defined  as  the  entanglement  of  each  pair,  divided  by  the  time  between  pairs. 
However,  such  a  naive  approach  fails  if  there  are  correlations  between  subsequent  pairs,  or  if  we  consider 
entangled  beams  of  radiation  that  cannot  be  naturally  decomposed  into  well-defined  pairs  of  particles.  In 
particular,  this  is  true  for  the  very  important  case  of  continuous  variable  (CV)  entangled  beams.  Although  many 
quantum  information  protocols  exploit  qubits,  with  their  discrete  state  space,  the  original  Einstein-Podolsky- 
Rosen  [1]  entanglement  involves  CVs,  and  CV  entanglement  has  many  modern  applications  [2-8]. 

In  this  article,  we  set  out  to  provide  a  general  definition  for  an  entanglement  rate  in  such  nontrivial 
situations.  It  will  turn  out  that  our  general  definition,  when  applied  to  stationary  Gaussian  CV  beams,  gives  rise 
to  a  frequency  integral  over  what  we  call  a  ‘spectral  density  of  entanglement’.  We  show  how  to  obtain  this  from 
the  two-point  time  correlators  of  the  entangled  beams,  using  a  suitable  additive  entanglement  measure  (the 
logarithmic  negativity  [9]).  For  one  of  the  most  common  situations,  we  also  provide  explicit  analytical 
expressions.  Our  definition  of  the  entanglement  rate  is  particularly  important  for  setups  where  the  output 
spectrum  is  arbitrary  (e.g.  containing  one  or  several  peaks).  This  is  a  widespread  case,  since  the  generation  of  CV- 
entangled  radiation  beams  is  often  enhanced  by  using  cavity  modes  (e.g.  [10, 1 1]  in  the  microwave  domain,  or 
[12, 13]  in  the  optical  domain).  Recently,  the  authors  of  the  first  experiment  on  spatially  separated  CV 
entanglement  in  superconducting  circuits  even  quantified  their  source  by  quoting  the  effective  number  of 
entangled  bits  per  second  [10],  estimated  from  the  bandwidth  of  the  circuit  and  the  entanglement  between  two 
modes.  Our  entanglement  rate  would  provide  a  precisely  defined  way  to  quantitatively  characterize  such 
situations. 

After  we  introduce  and  discuss  the  general  definition,  we  illustrate  our  entanglement  rate  by  applying  it  to  a 
four- mode  optomechanical  setup  that  allows  the  fully  resonant,  and  thereby  efficient,  generation  of 
entanglement.  The  rapidly  developing  field  of  cavity  optomechanics  focuses  on  the  dynamics  of  photons  and 
phonons  coupled  via  radiation  forces,  see  [14]  for  a  recent  review.  The  optomechanical  interaction  has  been 
predicted  to  produce  entanglement,  e.g.  between  optics  and  mechanics  [15-20]  or  between  light  modes  [21-31]. 
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Figure  1 .  (a)  A  source  of  two  entangled  beams,  where  the  aim  is  to  calculate  the  entanglement  between  beams  1  and  2 
interval  T.  (b)  Frequency  modes  that  are  correlated  (entangled)  in  the  limit  of  large  T. 

in  the  large  time 

Recently,  optomechanical  entanglement  was  demonstrated  experimentally  for  the  first  time,  in  a  microwave 
circuit  [32],  making  its  analysis  especially  timely. 

2.  The  entanglement  rate 

The  situation  we  have  in  mind  is  very  general:  a  source  emits  two  CV-entangled  beams,  described  by  bosonic 
fields  A\  ( t )  and  A2  (t),  see  figure  1  (a).  These  could  be,  for  example,  two  light  beams  of  different  polarization  or 
fields  propagating  along  two  different  waveguides.  Typical  sources  might  be  a  nonlinear  crystal  optical  resonator 
driven  by  a  pump  beam,  a  driven  optomechanical  cavity,  or  a  driven  on-chip  microwave  cavity  containing 
nonlinear  elements  like  Josephson  junctions. 

We  focus  on  the  important  regime  of  generating  stationary  Gaussian  CV- entanglement.  In  that  regime,  one 
treats  the  pump  as  classical  and  then  obtains  a  quadratic  Hamiltonian,  leading  to  Gaussian  statistics  of  the 
emitted  beams.  Because  of  the  pump,  that  Hamiltonian  will  be  time- dependent,  containing  terms  of  the  form 
dj  dj  e_KV  +  h.c.,  where  dc  would  be  a  cavity  mode  (our  analysis  also  applies  for  several  different  modes).  Here 
c Jp  is  the  pump  frequency  if  the  original  nonlinearity  was  of  the  x(2)  -type,  a]  a}  ap  +  h.c.,  whereas  uop  would  be 
twice  the  pump  frequency  for  a  x(3)  -type  nonlinearity  d }  a}  ap.  As  is  usual  in  such  situations,  it  will  be  most 
useful  to  switch  to  a  frame  rotating  at  the  frequency  uop /2,  such  that  the  Hamiltonian  becomes  time- independent 
and  we  are  dealing  with  a  stationary  problem.  In  that  new  frame,  u  =  0  relates  to  the  pump  frequency. 

Our  analysis  then  focusses  on  the  entanglement  properties  of  the  fields  Ay2  (t)  emitted  from  any  such  source. 
Being  Gaussian,  these  fields  are  completely  characterized  by  their  two-time  correlators.  The  details  of  the  source 
do  not  matter,  except  that  it  is  assumed  to  produce  Gaussian  beams  that  are  stationary,  i.e.  where  the  correlators 
only  depend  on  the  time-difference. 

At  this  point  we  note  that  in  some  situations  it  may  also  be  natural  to  consider  only  a  single  field  A  (t), 
propagating  along  a  single  waveguide.  Then,  frequency  components  centered  symmetrically  around  the  pump 
frequency  can  be  entangled,  and  they  may  afterwards  be  directed  to  two  different  output  ports  by  frequency¬ 
filtering.  For  such  a  situation,  we  can  still  apply  our  approach  if  we  define  A\  (t)  to  contain  the  positive  frequency 
(lj  >  0)  components  of  the  field  A  (t),  while  likewise  A2  (t)  would  contain  the  negative  frequency  (uj  <  0) 
components,  where  u  is  already  determined  with  respect  to  the  rotating  frame. 

Since  the  situation  is  stationary,  it  is  natural  to  try  and  define  an  entanglement  rate ,  i.e.  the  entanglement  per 
unit  time  emitted  from  the  source.  We  propose  to  do  this  in  the  most  natural  way  by  calculating  the  overall 
amount  of  entanglement  between  the  two  beams  in  a  long  time  interval  T  and  then  dividing  by  T,  in  the  end 
sending  T  to  infinity: 


r£  = 


lim 

T— oo 


E(T) 

T 


(1) 


This  definition  is  not  constrained  to  CV  entangled  beams  or  to  Gaussian  states.  It  only  requires  (a)  stationarity  of 
the  source,  (b)  an  entanglement  measure  that  is  additive  for  product  states,  and  (c)  a  finite  correlation  time  rc  for 
the  beams.  Given  such  a  finite  correlation  time,  the  fields  on  two  subsequent  time -intervals  of  length  T  rc  are 

not  correlated  to  a  very  good  approximation.  This  holds  because  even  though  there  maybe  remaining 
correlations  near  the  boundary  between  the  time-intervals,  these  are  appreciable  only  up  to  a  distance  of  order  rc 
from  the  boundary  and  can  be  neglected  in  the  limit  T  rc.  Due  to  the  additivity  (and  stationarity),  we  then  get 

E(2T)  ^  2E  (T),  such  that  the  entanglement  rate  calculated  for  time  intervals  of  sizes  Tor  2  T  is  the  same  (up  to 
the  small  corrections  which  we  can  neglect  in  the  limit  of  large  T). 
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From  now  on,  we  specialize  to  stationary  Gaussian  CV  entangled  beams.  We  will  use  the  logarithmic 
negativity  [4, 9]  ENas  an  entanglement  measure,  since  it  is  both  straightforwardly  evaluated  for  Gaussian  multi- 
mode  states  and  has  the  important  property  of  additivity. 

Let  us  now  consider  the  fields  As  ( t )  on  the  interval  [0,  T],  even  though  they  are  defined  for  all  times  t.  We 
can  define  discrete  frequency  modes  (s  =  1, 2): 

As{t)  =  f=  £  Awe-k.  (2) 

/V  >  ^  A  j 

The  normalization  is  chosen  such  that  the  ASyUJ  fulfill  bosonic  commutation  relations,  [ASyUJ,  Ay^> ]  =  6SyS>8UyU>, 
where  oo  =  j2ir/T  is  discrete,  with;  an  integer. 

We  now  want  to  calculate  the  full  logarithmic  negativity  EN  (T)  between  the  two  beams  on  that  time  interval, 
which  is  equivalent  to  the  entanglement  between  two  sets  of  harmonic  oscillators  [33].  In  our  case,  we  are  even 
considering  infinitely  many  harmonic  oscillators  As yU.  We  stress  that  the  entanglement  EN  (T)  is  (of  course) 
independent  of  our  choice  of  basis  for  each  of  the  beams,  as  a  different  choice  of  basis  amounts  to  implementing 
a  local  unitary  transformation.  The  correlations  between  the  two  beams  can  be  arbitrary,  except  that  they  are 
supposed  to  decay  beyond  some  correlation  time  rc  (which  is  true  in  any  reasonable  physical  situation).  As 
already  mentioned  above,  we  will  assume  that  T  rc,  such  that  the  correlations  between  subsequent  intervals 

of  size  T  can  be  neglected.  In  that  limit,  we  can  use  stationarity  to  show  that  only  the  following  types  of 
correlators  maybe  nonzero  in  the  present  situation,  up  to  corrections  that  are  small  in  rc/T:  [A] u  Ay  yU})  and 
{ASyUAyy-J)  (and  their  conjugates).  For  example,  we  find 

{^S,C0^s',C0')  ~  (As  Ay^—Qj ,  (3) 

where  we  have  defined  the  Fourier  transform  of  the  correlator:  (A5L  =  f +°°  dt  elu;t  (A  ( t)B  (0) ).  Likewise,  we 

J—  oo 

have  (ASyUJAyyUJ')  =  8^-^  ( AsAy  )UJ.  These  results  are  the  combined  outcome  of  stationarity  and  the  long  time 
interval  (i.e.,  T  ^  rc),  as  also  pointed  out  in  [34]. 

In  summary,  for  any  given  oo  >  0,  only  the  correlators  involving  the  modes  oo  and  —  oo  are  nonzero  (see 
figure  1(b)).  There  is  entanglement  EN  [oo]  between  the  two  modes  A1>u;  and  Ay_^  of  beam  1  with  the  two  modes 
A2 yCU  and  A2y-u;  of  beam  2.  Further  below,  we  will  show  how  EN  [oo]  can  be  calculated  from  the  correlators.  The 
overall  entanglement  En  (T)  between  the  two  beams  can  then  be  decomposed  into  a  sum,  because  of  the 
additivity  of  the  logarithmic  negativity: 

En(T)=J2En  M-  (4) 

u;>0 

We  note  that  EN  [oo]  in  the  sum  does  not  depend  on  T  (after  adopting  the  approximation  of  equation  (3)), 
although  the  number  of  summation  terms  does  scale  with  T  due  to  the  discretization  of  oo. 

The  sum  over  discrete  frequencies  can  be  converted  into  an  integral.  To  leading  order  at  large  T,  we  have 
Xw>0  EN  M  «  (27 r/T)  1  f  duj  En  [cifi.  Therefore,  we  can  identify  E]\j  [cid  as  the  spectral  density  of 
entanglement’4.  In  addition,  this  confirms  that  the  total  entanglement  indeed  grows  linearly  in  T for  sufficiently 
large  times.  As  a  consequence  we  finally  obtain  the  entanglement  rate: 

Ye  =  lim  =  f°°  ^En  M.  (5) 

oo  T  do  2ir 

This  shows  that  EN  [oo]  itself  maybe  interpreted  as  the  entanglement  rate  per  frequency  interval.  It  is  possible  to 
give  an  explicit  expression  for  EN  [oo]  in  terms  of  the  covariance  matrix  containing  the  correlators  [9].  The 
general  case,  involving  the  four  modes  A1>a;,  Ay_^  A2yUJ ,  and  A2j_ is  a  bit  cumbersome,  requiring  the 
symplectic  diagonalization  of  a  8  x  8  covariance  matrix.  However,  the  expressions  simplify  considerably  if  there 
is  purely  two-mode  squeezing,  i.e.  if  the  intra-mode  correlators  (A^A^^)  vanish.  In  that  case,  we  find 

En[oo]  =  E[oo]  +  E[-oo].  (6) 

For  positive  oo,  E  [oo]  is  the  entanglement  between  AljU;  and  A2j_ ^  while  E  [— oo]  is  the  entanglement  between 
Aiy_u  and  A2yCJ.  This  entanglement  between  two  opposite  frequency  components  has  been  analyzed  in  detail  in 
[34].  Especially  it  was  found  to  be  connected  with  the  composite  quadrature  variance  of  these  two  modes,  and 
therefore  could  be  measured  directly  in  experiments.  In  contrast  to  EN  [oof  the  density  E  [oo]  is  double-sided  (has 
contributions  both  at  negative  and  positive  frequencies).  We  note  that  in  the  following  we  will  also  commonly 
refer  to  E  [oo]  as  the  spectral  density  of  entanglement,  since  it  is  closely  related  to  EN  [oo].  Setting 

H+  =  (a!Ai,w>  +  \,n-  =  {aI_uA2,-J)  +  Y>andC  =  (A,uA-w}>  we  have: 

4 

Note  that  the  phrase  ‘entanglement  spectrum’  has  already  a  different  meaning,  which  is  why  we  prefer  to  use  ‘spectral  density  of 
entanglement’. 
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E[u]  =  max(0,  — ln(2?7_)), 

2r}_  =«+  +  «_-  7(«+  -  «-)2  +  4  l£l2  • 


(7) 

(8) 


As  an  aside  we  note  that  we  choose  to  work  with  the  natural  logarithm  in  our  discussion  (some  articles  use  log2, 
which  is  more  natural  for  discrete  qubits  [9]). 

In  this  special  case,  the  entanglement  rate  is  therefore: 


dcj 
27 r 


E  M. 


(9) 


3.  Relation  to  entanglement  between  wave  packets 

We  now  want  to  connect  our  general  result  to  previously  applied  approaches  for  quantifying  the  entanglement 
in  such  situations.  It  is  a  common  procedure  to  employ  normalized  mode  functions  (‘filter  functions’)  fs  ( t )  that 
have  the  shape  of  wave  packets,  in  order  to  define  two  harmonic  oscillator  modes,  one  for  each  beam: 

r+oo 

as=  /  d  tfs(t)As(t).  (10) 

J  —  OO 

Here  f +°°  \f  ( t )  |2d t  =  1.  The  entanglement  between  d\  and  a2  can  then  be  calculated,  e.g.  again  using  the 
logarithmic  negativity  as  an  entanglement  measure.  This  filtering  analysis  has  been  applied  to  several  settings 
[21, 26, 31],  for  a  detailed  explanation  see  e.g.  [21, 22].  The  catch  is  that  this  procedure  introduces  a  filtering  time 
r  (the  extent  of  the  wave  packets),  and  the  results  will  be  a  function  of  r,  which  is  usually  taken  to  be  arbitrary. 

How  would  one  loosely  define  an  entanglement  rate  based  on  this  procedure?  We  can  imagine  that  there  is  a 
stream  of  such  wave  packets,  with  a  spacing  of  about  r  (where  care  would  have  to  be  taken  to  define  them  to  be 
orthogonal).  A  simple,  though  phenomenological  approach  to  define  an  entanglement  rate  would  be  to  simply 
calculate  the  ratio  En/t. 

However,  it  is  clear  that  this  approach  is  not  systematic.  In  fact,  it  cannot  always  cover  the  full  entanglement, 
since  there  maybe  entanglement  in  components  of  the  beam  that  are  orthogonal  to  the  filter  functions  which  are 
employed.  In  addition,  there  is  some  arbitrariness  in  the  choice  of  filter  function  (and,  thus,  even  in  the 
definition  of  r).  Moreover,  one  may  have  situations  where  there  are  temporal  correlations  extending  beyond  r. 
Then,  the  entanglement  present  in  the  beams  maybe  underestimated.  If  one  tries  to  remedy  this  problem  by 
choosing  larger  r,  then  the  filter  bandwidth  shrinks  and  one  may  miss  entanglement  present  at  other 
frequencies. 

It  turns  out  that  an  approach  based  on  wave  packets  can  be  made  to  work,  but  only  if  one  constructs  a 
suitable  complete  basis  that  has  a  clear  physical  meaning.  It  is  more  systematic  than  the  naive  filtering  approach 
described  so  far  and  it  covers  all  the  entanglement.  In  addition,  it  can  be  related  to  a  direct  physical  prescription, 
and  we  will  see  that  it  leads  to  the  same  results  as  our  general,  basis-independent  definition  discussed  in  the 
previous  section. 

For  simplicity,  focus  on  the  situation  with  only  cross -correlations  (no  intra-beam  squeezing)  that  we 
discussed  at  the  end  of  section  2.  Imagine  one  sends  one  beam  through  a  frequency  filter  [a;,  u  +  8uo\  where 
8uo  =  27 r/r .  Likewise,  the  other  beam  will  be  sent  through  another  filter,  at  negative  frequencies 
[ — U,  —(jj  —  8uf\.  Now  construct  a  complete  set  of  orthogonal  wave  packet  modes  (‘Wannier  basis’)  with  a 
spacing  r  in  time,  which  are  able  to  fully  represent  the  filtered  beams  (see  figure  2).  As  we  show  in  the  appendix  A, 
the  logarithmic  negativity  En  between  two  such  wave  packet  modes  (one  in  each  beam,  at  equal  time-slots)  is 
related  to  E  [uj]  in  the  limit  r  rc: 

lim  E^j  =  E[cj].  (11) 

T— KX) 


As  a  consequence,  the  general  definition  of  the  previous  section  agrees  with  the  entanglement  rate  calculated 
from  such  a  wave  packet  picture.  This  wave  packet  approach  can  also  be  viewed  as  representing  the  following 
physical  procedure:  split  each  of  the  two  entangled  beams  into  many  frequency- filtered  output  beams,  where  the 
frequency  resolution  8u  =  2tt/t  has  been  chosen  fine  enough,  such  that  8lj  <C  1/tc.  The  rate  VE  quantifies  the 
total  entanglement  per  unit  time  contained  in  the  sum  of  those  streams  (since  it  was  defined  that  way  in  the 
previous  section).  Each  pair  of  wave  packets  (of  length  r)  can  in  principle  be  exploited  for  an  application  such  as 
CV  quantum  teleportation.  A  concrete  physical  measurement  of  the  entanglement  between  any  two  wave 
packets  could  be  performed  in  a  standard  way,  using  homodyne  measurements.  For  example,  the  local  oscillator 
can  provide  strong  pulses  that  are  shaped  in  the  form  of  the  wave  packet  modes  that  we  want  to  consider. 
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freq.  U) 


Figure  2.  Wave  packet  picture  that  can  be  used  in  understanding  the  entanglement  rate,  (a)  For  each  of  the  beams,  we  imagine  to 
frequency- filter  the  field,  with  a  filter  size  Suj.  (b)  In  time-space,  this  corresponds  to  wave  packets  spaced  by  r  =  In/Su;.  These  form  a 
complete  basis  on  a  time-frequency  grid.  The  use  of  a  wave  packet  basis  on  a  time-frequency  grid  in  the  context  of  quantum  noise  is 
reviewd  in  [35]. 


4.  A  specific  optomechanical  example,  and  its  implementation 


4.1.  Model 

We  illustrate  the  features  of  the  entanglement  rate  in  a  model  describing  the  effective  interaction  between  two 
localized  optical  modes  (d+  and  d_)  and  a  mechanical  mode  (b): 

H/H  =  —A  (d^d+  +  aid-)  +  Spb 

—  —  { (d+  +  al)b  +  (a^_  H-  d-)b},  (12) 

2 


Similar  Hamiltonians  have  been  studied  previously  in  the  context  of  entanglement  generation  between  light 
modes  [21-31].  For  instance,  Wang  and  Clerk  [28]  studied  intra- cavity  entanglement,  while  Tian  [29] 
investigated  also  the  stationary  output  entanglement.  We  note  that,  more  recently,  an  analysis  of  the  output 
entanglement  in  the  setup  with  A  =  6  =  0  (but  with  unequal  couplings  for  up-  and  down -conversion)  was 
provided  in  [30].  The  authors  gave  analytical  insights  into  the  spectral  density  of  entanglement  E  [cj].  It  was 
found  that  equal  couplings  can  lead  to  a  much  higher  maximum  value  of  E  [uj],  which  shows  up  at  a;  =  0.  But 
the  entanglement  centered  around  u  =  0  drops  rapidly  when  some  finite  frequency  bandwidth  is  considered. 

To  relax  the  sensentivity  to  bandwidth,  they  introduce  an  optimal  delay  time  between  the  two  filtered  modes. 
Here,  we  only  discuss  zero-bandwidth  entanglement.  To  optimize  the  output  entanglement,  we  therefore  choose 
equal  couplings  for  the  up-  and  down- scattering. 

Before  proceeding,  we  note  how  to  implement  this  model  using  three  equidistant  optical  modes,  enhancing 
the  efficiency  beyond  previous  suggestions.  The  optical  mode  spacing  /  is  nearly  resonant  with  the  vibration 
frequency  El,  with  a  frequency  mismatch  8  =  El  —  J  (figure  3(a)).  A  laser  drives  the  center  optical  mode  at 
with  a  detuning  A  =  uL  —  a;0.  An  optomechanical  interaction  of  the  kind  %0d^d0(^  +  P),  and  similarly  for 
d_,  scatters  photons  up  and  down,  into  modes  d+  and  d_,  while  simultaneously  destroying  (creating)  phonons. 
When  a  phonon  is  virtually  emitted  and  re-absorbed,  an  effective  four-wave  mixing  process  is  induced, 
generating  a  pair  of  a+  and  d_  photons  out  of  a  pair  of  d0  photons.  Thus,  two-mode  vacuum  squeezing  (EPR 
entanglement)  is  produced.  We  assume  that  the  drive  is  strong  and  we  can  replace  d0  by  the  coherent  amplitude 
a  =  (d0).  This  yields  the  Hamiltonian  (12)  with  g  =  2 gQ  a ,  provided  we  choose  frames  rotating  at  uoL  ±  /  for 
the  modes  d± ,  and  at  /for  b .  Moreover,  only  nearly  resonant  terms  are  kept,  which  is  allowed  if  /,  El  g,  k, 
where  k  is  the  optical  intensity  decay  rate. 

Possible  experimental  implementations  include  a  membrane-in-the-middle  setup  tuned  to  a  point  with 
three  equidistant  modes  [36-38]  (figure  3(b))  or  coupled  optomechanical  cells,  e.g.  in  an  optomechanical  crystal 
[39, 40]  (figure  3(c)),  see  the  appendix.  Such  a  triply  resonant  setup  enhances  the  efficiency:  as  compared  to 
previous  suggestions  with  only  one  resonantly  driven  mode,  generating  entangled  Stokes  and  anti- Stokes 
sidebands  (similar  to  [21]),  one  wins  a  factor  ElA  /  (d2  +  (ft/2)2)2  ^  1  in  the  intensity  of  the  entangled  output 
beams,  while  compared  to  setups  with  two  optical  ouput  modes  (e.g.  [28, 29]),  one  wins  a  factor  (2El / k)a  1, 

for  fixed  input  laser  power  and  A  =  0. 

We  use  standard  input-output  theory[41]  for  our  analysis: 


aj  —  [H,  ad  aj  ^fiEay m(t), 

a  2 

£  =  T[A,  b]  -  jb  -  A/r^t). 


(13) 
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Figure  3.  (a)  Level  diagram:  the  optical  level  distance  almost  matches  the  vibrational  frequency  so  that  phonon-mediated  transitions 
between  the  optical  levels  occur,  (b)  Schematic  for  a  possible  implementation  in  a  membrane-in-the-middle  setup;  three  of  the  optical 
modes  constitute  an  equidistant  triplet  for  a  suitable  membrane  position  and  tilt  angle,  (c)  Schematic  for  the  possible  implementation 
in  an  optomechanical  crystal,  with  one  incoupling  and  one  outcoupling  waveguide. 


Here  T  is  the  mechanical  damping  rate  and  j  =  =L  As  usual,  (b^n  ( t)bm  (0))  =  n^S  ( t )  and  (bin  (t)b^n  (0))  = 
(n^  +  1)<5(0,  with  n^  =  (exp  (Ml/{k^T))  —  l)-1  the  thermal  occupation,  and  likewise  for  n  (but 
without  thermal  noise).  Solving  equation  (13)  and  employing  dyout  (t)  =  dj^n(t)  +  Vftd;(0>  we  find  the 
linear  relation  between  output  and  input  fields  in  terms  of  a  scattering  matrix. 

Coming  back  to  our  general  definition  of  the  entanglement  rate,  we  would  consider  d+jOUt  ( t )  as  the  first 
beam  A\  (t)  and  d_>out  (t)  as  the  second  beam  A2  ( t ).  It  is  not  difficult  to  show  (and  can  be  confirmed  by  direct 
calculation)  that  there  is  no  intra-beam  squeezing,  i.e.  (Ai^Ai^J)  =  0  and  likewise  for  beam  2.  Thus,  we  want 
to  employ  the  formulas  equations  (7)  and  (9)  in  order  to  find  the  contributions  to  the  spectral  density  of 
entanglement  and  the  entanglement  rate. 

To  calculate  E  [a;],  we  need  the  correlators  of  the  two  beams.  Entangled  photon  pairs  are  emitted  at  physical 
frequencies  ljl  ±  J  ±  l u,  corresponding  to  db  uj  in  our  rotating  frame.  We  have  to  evaluate  correlators  like 
j  °°  eluJt  (d+,out  (t)d_jOUt  (0) )  d t  as  shown  in  appendix  E. 


4.2.  Results 

We  first  address  some  general  features.  The  system  can  become  unstable  (figure  4(a)),  both  towards  mechanical 
and  optical  oscillations.  The  optical  stability  boundary  is  approximately  given  by  (A  +  g2/26)A  +  (ft/2)2  =  0. 
Eliminating  the  mechanical  mode  (for  6  ft,  T,  g,  |  A  |),  we  obtain  an  effective  optical  model,  which  yields  the 

photon  pair  creation  rate  rpairs=  y  (ft/2)2  +  ^  A  +  A  j  .  This  diverges  at  the  optical  stability 
boundary. 

We  now  focus  on  the  spectral  density  of  entanglement  E  [cj]  that  characterizes  the  output  beams,  and 
especially  the  entanglement  rate.  In  contrast  to  the  intra-cavity  entanglement  (discussed  in  [28, 29]),  we  find  that 
E  [uj]  is  not  bounded.  This  is  similar  to  the  difference  between  intra-cavity  and  output  squeezing  [42-44] . 
Numerical  plots  of  E  [cj]  for  the  Hamiltonian  (12)  have  been  shown  so  far  [29]  only  for  the  special  case 
A  =  6  =  0  and  asymmetric  optomechanical  couplings.  Entanglement  of  temporal  modes  was  also  discussed 
for  a  pulsed  scheme  [3 1]  in  the  case  6^0. 

The  output  entanglement  grows  significantly  at  the  stability  boundary  (figure  4(b)),  even  diverging  in  the 
effective  optical  model.  This  is  typical  near  an  instability  but  it  comes  at  the  price  of  linewidth  narrowing, 
reducing  the  entanglement  rate. 

In  order  to  appreciate  this,  we  now  discuss  the  output  intensity  spectrum,  S+  (u)  =  J e~luJt 

(d^ out  (0  <2+, out  (0) )  df ,  in  figures  5(a)-(c).  We  expect  that  any  mechanical  noise  contributing  to  the  optical 
output  is  deleterious  for  entanglement,  which  is  confirmed  by  explicit  calculation.  Thus,  we  plot  the  spectrum  in 
an  instructive  way,  distinguishing  optical  and  mechanical  contributions,  as  obtained  from  the  linear  relation 
d+>out(a;)  =  ...fl+jin(o;)  +  ...(a/ in)(o;)  +  (a;).  There  are  typically  two  peaks,  separated  by  8  =  El  —  J 
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Figure4.  (a)  Stability  diagram  of  the  model,  versus  frequency  mismatch  <5  and  laser  detuning  A,  for  g/n  =  5  and  T/k  =  10-3. 
Stability  boundary  in  the  effective  optical  model:  red  dashed  line.  Parameter  values  corresponding  to  the  (stable)  blue  points  will  be 
studied  below,  (b)  Diverging  entanglement  from  the  effective  optical  model  at  the  boundary  of  stability. 
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Figure  5.  The  output  spectrum  S+  (c j)  for  the  beam  from  the  upper  optical  mode,  plotted  for  three  typical  situations  (a)-(c)  at  a  fixed 
temperature  («t h  =  50),  with  g/K  =  5  and  T/k  =  10-3.  Blue  versus  brown  distinguish  optical  versus  mechanical  noise 
contributions,  (a)  resonant  drive,  off-resonant  mechanical  mode  (A  =  0,  6  =  10k);  (b)  near  the  optical  instability 
(A  =  —0.2 k,  6  =  10k);  (c)  doubly  resonant  (A  =  6  =  0).  Note  the  differences  in  vertical  and  horizontal  scales  for  the  peaks,  (d)-(f) 
The  spectral  density  of  entanglement  E  [cj],  corresponding  to  (a)-(c).  Additional  curves  in  (f)  correspond  to  T/k  =  5  x  10-2 
(dashed)  and  T/k  =  10-3  (solid),  for  n±  =  0  (orange)  and  n±  =  50  (black),  at  fixed  cooperativity  C  =  £2/kT  =  2.5  x  104. 


and  containing  primarily  optical  or  mechanical  noise,  respectively.  Near  the  optical  instability  (figure  5(b)),  the 
optical  peak  gets  strong  and  narrow. 

When  the  optical  mode  spacing  matches  the  mechanical  frequency  (6  =  0)  and  the  laser  is  on  resonance 
(A  =  0),  the  two  peaks  merge  and  have  a  narrow  linewidth  set  by  the  mechanical  damping  rate  T  (figure  5(c)).  It 
will  turn  out  that  the  entanglement  rate  is  maximized  near  this  point.  This  is  entirely  counterintuitive:  One 
might  expect  that  at  6  =  0  mechanical  noise  is  injected  into  the  output  beams,  destroying  entanglement. 
However,  we  find  that  the  optical  noise  can  completely  overwhelm  the  mechanical  noise  for  strong  driving, 
when  the  cooperativity  is  sufficiently  large,  C  =  g2 / k T  nt h. 

The  spectral  density  E  [uj]  is  shown  in  figures  5(d)— (f).  Typically,  E  [cd]  is  maximal  at  the  optical  peak  near 
u  =  0.  For  the  important  case  A  =  8  =  0,  we  find  an  analytical  expression,  E  [cj  =  0]  =  —  In  (2rj_),  where 


El-  —  4  C^C  +  nth  +  —  j 


-2CJ1  +  4  C 


(c+»th  +  i) 


(14) 


depends  on  both  driving  (via  C)  and  temperature.  We  checked  that  choosing  different  optomechanical 
couplings  for  the  two  modes  d±  will  not  increase  E[u  =  0],  in  contrast  to  the  intra-cavity  case  [28]. 
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Figure  6.  (a)  Spectral  density  of  entanglement,  maximized  over  frequency,  maxw£  fij],  and  (b)  total  entanglement  rate 
T e  =  Je  [cj]  du;  j  2tv  (in  units  of  «),  both  shown  as  a  function  of  frequency  mismatch  6  and  laser  detuning  A,  for  n±  =  0,  g/n  =  5 
and  T/k  =  10-3.  (c)  The  full-width-half-maximum  (FWHM)  of  the  peak  in  E  [uj]  as  a  function  of  T/ac,  for  fixed  cooperativity 
C  —  2.5  x  104  (solid)  and  C  =  103  (dashed),  at  6  =  A  =  0.(Blue:  n±  =  0,red:  nth  =  50).  (d)  Temperature  dependence  of  Te 
(6  =  A  =  0  and  F/k  =  10-3.) 


In  figures  6(a)  and  (b),  we  compare  the  maximum  of  the  spectral  density  of  entanglement,  £max  =  maxw£  [a;], 
and  the  entanglement  rate  Te  =  Je  [uj]  dujflir.  While  £max  becomes  large  near  the  optical  boundary  of  stability, 
the  entanglement  rate  there  remains  small,  due  to  the  narrow  bandwidth.  This  will  be  a  general  feature  in  many 
similar  systems.  Rather,  VE  is  optimal  near  6  =  A  =  0.  We  found  that,  as  the  mechanical  damping  rate  increases, 
the  optimum  shifts  away  from  A  =  6  =  0. 

The  entanglement  rate  depends  on  the  full  shape  of  E  [uj],  in  particular  the  peak  width(s).  Crucially,  those 
widths  are  distinct  from  those  in  the  output  spectrum,  due  to  the  nonlinear  (logarithmic)  dependence  ofE  on 
parameters.  For  example,  in  the  case  A  =  6  =  0  of  greatest  interest,  the  peak  width  (see  figures  5(f)  and  6(c))  is 
not  set  by  the  small  mechanical  linewidth  T,  unlike  for  the  output  spectrum  discussed  above.  Thus,  the  values  of 
the  entanglement  rate  VE  for  the  parameters  explored  here  are  much  larger,  of  the  order  k  •  ma x^E  [a;].  For 
increasing  temperatures,  VE  decreases,  but  only  slowly,  indicating  robust  entanglement  in  this  setup:  VE  oc  n^1, 
with  the  prefactor  set  by  the  cooperativity  C  (figure  6(d)). 


5.  Conclusions 

We  have  introduced  an  entanglement  rate  as  a  quantitative  measure  for  the  CV  entanglement  production  per 
unit  time  in  setups  involving  resonant  modes.  The  definition  is  natural,  in  that  it  simply  characterizes  the  total 
entanglement  between  two  beams  within  a  time-interval  of  size  T,  in  the  limit  T  — >  oc.  In  principle,  it  is  also 
more  general  than  the  Gaussian  CV  case  studied  in  the  present  manuscript. 

Moreover,  we  have  studied  an  optomechanical  setup  with  one  mechanical  and  three  optical  modes  that 
allows  fully  resonant  production  of  optical  entanglement.  The  spectral  density  of  entanglement  and  the  overall 
entanglement  rate  are  optimized  for  different  parameter  choices.  The  concept  introduced  here  should  be  useful 
for  analyzing  setups  in  other  domains,  like  cavities  with  a  Kerr  medium  [12, 13]  or  microwave  resonators  with 
nonlinearities  [10, 1 1]. 
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Appendix  A.  Wave  packet  basis:  general  scheme 

In  the  main  text,  we  claimed  that  our  general  entanglement  rate  can  be  connected  to  an  approach  of  frequency- 
filtering  and  decomposing  the  filtered  beams  into  wave  packets.  The  natural  way  to  make  these  notions  precise  is 
by  using  the  Wannier-basis  of  wave  packets  fm  n  ( t )  that  live  on  a  regular  grid  both  in  time,  tn  =  nr ,  and  in 
frequency- space,  uom  =  mduo,  where  duo  =  27t/t.  The  Fourier-transform  Fm,n  ( uo )  =  J n  (t)dt /  J2tt  of 
these  wave  packets  is  nonzero  only  in  the  interval  a;  E  ^(m  —  (m  +  where 

Fm,n  (<^0  =  exp  (ic otn)/ > Jduo .  The  Fm^n  ( uo ),  and  therefore  also  the  fm  n  (t),  form  an  orthonormal  basis.  That  makes 
it  possible  to  decompose  uniquely  the  output  field  into  modes  defined  by  these  wave  packets: 

As(t)  =  where  a$n  =  ff*  n(t)As(t)  dt  annihilates  a  photon  in  mode  (m,  n),  and 

[a$n,  dj--]  =  In  atypical  application  with  a  stationary-state  source  producing  beams  via 

parametric  down- conversion  or  four- wave  mixing,  energy  conservation  dictates  that  pairs  of  frequencies  uo\  and 
uo2  are  entangled  only  when  they  obey  a;totai  =  uo\  +  uo2  (in  the  limit  r  — >  oo).  In  the  rotating  frame  adopted  in 
the  main  text,  uotot3L \  would  be  zero. 

We  now  calculate  the  logarithmic  negativity  En  between  two  wave  packets  of  the  two  output  fields,  at  the 
same  time  slot  tn  and  in  suitable  frequency  slots  uo  and  —uo.  For  simplicity,  we  will  assume  a  situation  with  purely 
inter-beam  entanglement  (no  intra-beam  squeezing),  as  discussed  at  the  end  of  section  2.  A  similar  construction 
would  apply  to  the  more  general  case,  even  though  then  it  would  become  necessary  to  consider  the  entanglement 
between  two  wave  packets  at  uo  and  —uo  of  beam  1  with  two  wave  packets  at  uo  and  —uo  of  beam  2. 

To  prepare  for  our  definition  of  the  entanglement  rate,  we  have  to  discuss  the  dependence  of  E Jj  on  the  filter 
time  r,  especially  for  the  limit  r  rc,  where  r  is  much  longer  than  the  physical  correlation  time  rc  of  the  source 

( r~l  is  the  width  of  the  spectral  peaks).  First,  this  ensures  that  there  will  be  no  correlations  between  wave  packets 
located  at  different  time  slots,  such  that  we  capture  the  full  amount  of  entanglement.  Second,  we  will  now 
explain  in  this  wave  packet  picture  why  En  tends  to  a  well-defined  limit  for  r  — >  oo,  which  is  consistent  with 
direct  calculations  [21]. 

Consider  enlarging  the  filter  time  r  by  a  factor  M,  shrinking  the  filter  frequency  interval  by  duo'  =  duo /M.  In 
effect,  the  new  wave  packets  of  size  r'  =  Mr  encompass  M  of  the  old  wave  packets.  That  this  coarse- graining 
keeps  the  correlators,  and  thus  the  entanglement  Etn,  unchanged  can  be  understood  already  from  a  very  simple 
consideration.  Take  a  suitably  normalized  sum  of  M  operators,  X  =  Mrxd2Y^=  \  Xn  and  likewise 
Y  =  M~l/2Y^=i  Yn-  Then  the  correlator  between  these  ‘averaged5  operators  will  equal  the  original  correlator: 
(XY)  =  (Xn  Yn).  This  holds  provided  there  have  been  no  cross- correlations  and  (Xn  Yn)  is  independent  of  n.  The 
same  logic  applies  to  our  case,  where  the  new  modes  are  properly  normalized  averages  over  the  old  modes: 
am^nt  =  ^2nKmf(n  —  n'M)  dm>„,  with  ^encoding  the  overlap  between  the  two  basis  sets.  The  detailed  structure 
of  Kis  not  important  for  our  argument,  but  essentially  K  is  nonzero  only  in  a  range  of  size  M,  for 
|  n  —  n'M  |  <  M,  and  it  is  normalized:  | Km>(n  —  n'M)  \2  =  l.As  shown  in  the  next  section,  this  ensures  a 

well-defined  limit  E^°°. 

The  entanglement  rate  for  a  given  frequency  slot  m  may  now  be  defined  naturally  as  E^ >m/r,  the  ratio 
between  the  entanglement  ETN,  m  contained  in  a  pair  of  wave  packets  at  that  frequency  and  the  time  r  between 
those  wave  packets.  (Here  m  would  denote  the  index  for  beam  1  to  which  corresponds  uniquely  an  index  m2  for 
beam  2,  as  explained  above.)  As  the  logarithmic  negativity  is  additive,  it  makes  sense  to  add  up  Ek  m  for  all  the 
small  frequency  intervals  into  which  the  emission  of  the  source  has  been  decomposed. 

We  have  checked  by  direct  calculation  that  limr^00  defined  here  co-incides  with  the  E  [uo]  defined  in 

the  main  text  (where  uo  =  mduo  is  kept  fixed  in  the  limit  r  — >  oo,  by  adjusting  m).  This  is  despite  the  fact  that 
EL  m  is  defined  with  respect  to  Wannier  basis  modes,  while  E  [uo]  was  defined  from  the  Fourier  components  of 
the  field  on  a  finite  time-interval  of  length  T  rc. 
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We  thus  have  arrived  again  (in  a  different  route),  at  the  total  entanglement  rate: 


I'e  =  lim 


~'N,m 


f+°° 

’Joe  2i r 


(15) 


We  have  used  8w  =  2r /t  to  convert  the  sum  into  an  integral.  In  the  limit  r  — >  oo,  the  details  of  the  Wannier 
basis  have  become  unimportant.  E  [a;]  can  now  be  calculated  using  any  filter  function. 


Appendix  B.  More  details  on  the  Wannier  basis 


In  this  section,  we  present  some  more  technical  details  on  our  wave  packet  based  discussion  of  the  entanglement 
rate.  As  defined  in  the  main  text  and  mentioned  in  the  preceding  section,  we  start  from  a  complete  orthonormal 
basis  of  packets  that  are  localized  both  on  a  time-grid  (tn  =  nr)  and  a  frequency- grid  (ujm  =  mfiuj),  where 
r  =  hr/Sco.  In  frequency  space,  these  basis  functions  are  defined  as 


Fm,n  (^0 


(16) 


if(m  —  1/2 )6uj  ^  lu  <  (m  +  1/2)  fe,  and  Fm>n(c<;)  =  0  otherwise.  In  time-space,  this  corresponds  to  the 
well-known  Wannier-basis  of  sine-shaped  wave  packets 

=  f  Fm,n(u)e-iuJt-E=  =  fm(t  ~  tn), 

where 


/m« 


1  c-ia^sin(&J*'/2) 

Vr  8ut/2 


(17) 


Upon  temporal  coarse- graining,  we  combine  M  old  wave  packets  into  one  new  one,  and  at  the  same  time  the 
frequency-resolution  becomes  refined:  r'  =  Mr  and  8uo'  =  8uj/M.  Thus,  the  new  frequency  index  m'  can  be 
thought  of  as  a  combination  of  the  old  index  m  and  another  index  l  =  0  ...  N  —  1  that  splits  the  old  frequency 
interval  into  Mpieces.  The  interval  belonging  to  m'  is  thus:  u  E  [(m  —  1/2 )8u  +  18cj',  (m  —  1/2 )8u  + 

(I  +  1)  &/[.  The  definition  of  F'mi  ,  (u)  on  this  interval  reads  like  the  old  one,  except  for  the  obvious 
replacements:  fv  n/  (u)  =  eluJtn  f  V 8uj'  ,  with  t'nt  =  n'r'  =  n'Mr.  Both  the  old  and  the  new  basis  are  complete. 

We  now  want  to  obtain  the  overlap  integrals  that  relate  the  old  basis  to  the  new  one.  We  find: 

(m’,  n’\m ,  n)  =  J F'*m'y(cv)Fm?„(u) Au  =  Km:(n  -  n'M ),  (18) 


with 


Km'(k ) 


k  =  0: 
k  ^  0: 


i 

4m 


Vm 


<Pmk  -  1 
2irik 


(19) 


Note  that  the  overlap  K  obviously  depends  on  the  frequency  index  m'  only  via  the  refinement  index  /  in 
m'  =  (m,  /),  and  that  K  =  0  if  we  were  to  calculate  the  overlap  for  any  m'  that  is  not  part  of  the  original 
frequency  interval  defined  by  m.  The  shape  of  the  overlap  as  a  function  of  the  temporal  distance  n  —  n'M  is  that 
of  a  sine  function  with  a  decay  scale  set  by  M.  One  can  confirm  that  the  overlap  matrix  elements  calculated  here 
are  normalized,  \Km'(n  ~  n'M)\2  =  1. 

When  we  think  of  the  situation  with  two  entangled  beams,  we  imagine  each  of  them  is  defined  by  its  own 
fluctuating  output  field,  Aa  (t ),  where  a  =  1,2  denotes  the  beam.  Each  of  those  can  be  decomposed  into  the 
kind  of  Wannier  basis  defined  here,  and  we  choose  the  same  filter  time  r  for  each  of  them.  Regarding  the 
frequencies,  we  want  to  exploit  the  fact  that  in  a  typical  steady- state  situation  (like  in  parametric  down- 
conversion  or  four-wave  mixing),  it  is  pairs  of  frequencies  that  are  entangled.  Thus,  to  each  frequency  of 
beam  1  belongs  an  entangled  frequency  cu2  of  beam  2  (with  a;tota  i  =  uj\  +  o;2).  To  simplify  the  subsequent 
notation,  we  want  to  shift  and  revert  the  frequency  scale  of  beam  2  such  that  the  two  mutually  entangled 
frequencies  are  always  both  denoted  by  the  index  m.  In  other  words,  while  uj\  =  mScj,  we  have 
<^2  =  cc;totai  —  mbuj.  We  note  that,  after  going  into  a  rotating  frame  that  makes  the  Hamiltonian  time- 
independent  (as  in  the  main  text),  we  obtain  a;tota i  =  0.  In  addition,  it  turns  out  that  due  to  this  matching 
between  opposite  frequencies,  the  basis  functions  for  beam  2  have  to  be  changed  accordingly,  and  the  basis 
transformation  for  beam  2  is  effected  by  K *  instead  of  K. 

We  now  want  to  confirm  (as  indicated  in  the  main  text),  that  such  a  coarse -graining  does  not  change  the 
entanglement  E  between  wave  packets,  provided  we  are  already  in  the  regime  of  sufficiently  large  filter  time, 
r  rc.  To  this  end,  we  just  have  to  show  that  the  correlators  between  modes  do  not  change  upon  coarse- 
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graining,  i.e.  we  want  to  show  in  particular 


\um,num,n/  \um' ,n'um' ,n' /  > 


(20) 


where  the  modes  on  the  right-hand  side  are  the  temporally  coarse-grained  ones,  and  ujm'  lies  within  the  interval 
defined  by  m.  The  precise  location  of  ujm'  within  that  interval  will  not  matter,  since  6uj  <C  l/rc,  so  the  spectrum 
of  the  source  is  already  flat  on  that  scale.  In  addition,  we  note  that,  for  the  steady-state  situation  we  assume  here, 
neither  the  left-hand  side  nor  the  right-hand  side  of  (20)  actually  depend  on  the  time-point  n  or  rc',  respectively. 

Employing  the  overlap  calculated  above,  we  find: 


/a'(D  ^ '(2)  \ 

\am\n,am\n’l 

=  XX'Cm  -  n'M)K*,(n2  -  n'M)(&2?„A) 

nhn2 

=  £ \K„'  (n  -  n'M)  |2  (O Z )  =  (0®«)  ■  CD 

n 

Here  we  have  used  that  a^Pn,  =  X„,-^m'(«i  -  n'M)  a ^  and  a'ffn,  =  T,niKm'(n2  ~  In  going  to 

the  second  line,  we  exploited  the  fact  that  different  time  slots  are  already  uncorrelated  (r  rc).  We  also  used  the 

normalization  of  Kin  the  last  step. 

Regarding  other  correlators,  such  as  we  can  say  the  following:  for  the  model  discussed  in  the 

main  text,  they  can  be  confirmed  to  be  zero  by  explicit  calculation.  For  more  general  models,  where  however  the 
Hamiltonian  can  still  be  brought  to  a  time-independent  form  by  a  proper  choice  of  rotating  frame,  we  find  that 
stationarity  dictates  that  (d(1)  (uj)  (d(2)  (u/) )  oc  6  (uj  +  uj'),  and  likewise  for  all  other  correlators.  In  evaluating  a 
correlator  of  amplitudes  in  the  Wannier-basis,  like  we  are  effectively  looking  at  correlators  of  the 

type  (d(1)  (o;)(d(2)  (a/)/),  with  uj'  ~  —  uj  [by  our  definition  given  above,  m  enters  the  frequency  uj2  with  an 
opposite  sign].  Since  (d(2)  (a;') 'f  =  (a(2)1")(— a/),  that  correlator  equals  (d(1)  (uj) (d^2^)(— ur))  oc  6(u  —  uor).  By 
virtue  of  a/  w  —uj,  this  is  zero  automatically.  (Note  that  formally  we  have  to  exclude  the  single  point  u  =  0  in 
this  argument,  which  would  physically  equal  the  incoming  laser  frequency  in  the  case  of  four- wave  mixing,  and 
the  corresponding  small  frequency  slot.)  For  the  correlators  involving  quantities  of  the  same  beam,  it  is  rather 
correlators  of  the  type  (d^rAm^n)  that  are  nonzero,  while  correlators  like  (a^na^n)  are  zero,  and  the  proof  is  in 
analogy  to  what  was  shown  above.  In  summary,  the  entanglement  E  will  go  to  a  well-defined  limit  as  r  — >  oo 
(t  rc).  (As  mentioned  above,  in  taking  this  limit,  we  assume  the  frequency  index  m  is  re-adjusted  such  that 

the  center  frequency  of  the  corresponding  slot  is  held  fixed  at  some  given  u,  thus  arriving  at  E  [uj]  in  the  limit 
t  — >  oo.)  The  actual  calculation  of  E  [uj]  is  discussed  below. 


Appendix  C.  Hamiltonian  and  implementation  with  three  coupled  optomechanical  cells 


In  this  section,  we  give  a  derivation  of  the  Hamiltonian  (2)  and  discuss  its  implementation.  We  consider  three 
equidistant  optical  modes  with  a  splitting /coupled  to  the  same  mechanical  mode  b  with  frequency  Q  via 
radiation  pressure.  One  of  the  optical  modes  (here  called  a0)  is  driven  by  an  external  laser  at  frequency  04.  Such  a 
setup  in  general  can  be  described  by  the  following  Hamiltonian  (h  =  1),  where  for  brevity  we  do  not  display  the 
coupling  to  the  dissipative  environment: 


H=  ^2  u}qa]aq  +  Qb*b  -  g^qa^,aq(b  +  P)  +  A(a0ei(oJLt+,l>)  +  h.c.).  (22) 

q—zL,0  q,qr=Jc,  0 

Herecc;±  =  uj0  ±  /,  and  represents  the  generic  (Hermitian)  optomechanical  coupling  matrix.  We  now 

q  ,q 

assume  that  the  frequency  mismatch  6  =  Q  —  /  is  much  smaller  than  the  mechanical  frequency,  i.e.,  |<5|  H 
and  that  the  mechanical  sidebands  are  resolved,  i.e.,  k  <C  H,  where  k  is  the  optical  damping  rate.  After 
transforming  to  the  rotating  frame  with  respect  to  H0  =  (ujl  +  /)a^d+  +  (wL  —  J)ala _  +  ujLa^do  +  Jb^b, 
and  neglecting  rapidly  oscillating  terms  rotating  at  d=/  by  a  rotating  wave  approximation  (RWA),  we  find 


Hr 


■RWA 


E 

q—^z,0 


A  a\aq 


+ 


bb'b 


{(^)°o«+«o  +  ga]Ala-)b  +  h.c.J  +  A (a0e^  +  a£e  i<j>) 


(23) 


where  A  =  wL  —  uj0  is  the  laser  drive  detuning.  For  a  sufficiently  strong  laser  drive,  we  can  linearize  the 
dynamics  by  replacing  d0  by  a  complex  number  a  =  —  - .  Thus  we  find  the  Hamiltonian  (2)  as  given  in 

the  main  text 


H\in  —  — A  (fT^$_|_  T  cE_  uE)  T  6b  b  —  T  U—)b  T  h.c. },  (2 4) 

provided  that  the  couplings  to  the  £+’  and c— ’  mode  turn  out  to  be  equal,  i.e.  g/2  =  =  g, Without 

loss  of  generality,  we  assume  that  g  is  real- valued.  The  RWA  made  above  is  valid  when  \g\  <C  fh  It  maybe 
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unavoidable  that  there  is  a  small  asymmetry  in  the  two  optomechanical  coupling  strengths,  but  since  the 
entanglement  generation  involves  both  transitions  simultaneously  this  does  not  make  a  crucial  difference  (we 
have  confirmed  this  for  the  output  entanglement  which  we  discuss).  A  more  detailed  account  on  some  further 
aspects  of  asymmetric  optomechanical  coupling  strengths  in  the  context  of  this  kind  of  Hamiltonian  maybe 
found  in  [28, 29],  where  it  is  pointed  out  that  for  the  case  of  intra-cavity  entanglement  (which  is  not  our  concern 
here)  having  such  an  asymmetry  can  actually  be  beneficial. 

We  now  turn  to  deriving  the  Hamiltonian  (22)  for  a  concrete  setup  consisting  of  three  coupled 
optomechanical  cells.  In  each  of  the  cells,  we  assume  a  standard  (local)  coupling  between  photons  and  phonons. 
The  microscopic  Hamiltonian  reads  (h  =  1) 


hfcell 


^  ^  t  ^  ^  t  /V 

^{CJ0 fl/flz  +  fib/  h  -  g0*i*i(bi  +  bi)} 

i=i 


(K\ a i  a2  +  +  h.c.), 


(25) 


where  the  index  /  =  1 ,  2,  3  runs  over  the  three  sites.  The  second  term  describes  photon  tunneling  between 
different  sites.  It  can  be  taken  into  account  by  introducing  optical  normal  modes,  as  defined  by 

a±  =  (  K]U]  +  =  =F  a2 1  /  V2 ,  d0  =  ~  KlU3  with  eigenfrequencies  uo±  =  cj0  ±  JKi  +  fC22,  respectively.  In 

+  )/  +  v 
terms  of  the  normal  modes,  the  Hamiltonian  can  be  written  as 

3  -t-  -t 

HCell  =  E  +  YNhl  6  “  f  («-  “  “+f  («-  ~  «+)(&2  +  fa) 

q=±,0  1=1 

~  +  K\d-  +  42K2dof  {K\d+  +  K\d-  +  j2K2do)(bl  +  £1) 

2  (A]  +  a2  ; 

—  r2)  +  K2d—  —  {K2a+  +  K2a _  —  V2X1a0)(b3  +  £3).  (26) 

2  (A]  +  A2  ) 


This  Hamiltonian  takes  essentially  the  same  form  as  the  Hamiltonian  (22)  given  above.  We  assume  that 
*  =  n  -  Jk[Tk2  <  £2  and  ft  <C  ^2,  transform  to  a  rotating  frame  with 

=  E^±,oH  -  +  E?=  1  V^i2  +  ^2  bi  bh  and  apply  a  RWA  to  find 


Heff  —  ^  ^ 


<2  a<2 


<2=i,0 


E^/6 

/=1 


KiK2 


Kf  +  Ki 


b 

(u^_Uq  +  U—Uq  )  — 


V2 


h.c. 


(27) 


After  adding  an  external  laser  drive  for  the  d0  mode,  moving  into  a  frame  rotating  with  the  drive  frequency  and 
applying  standard  linearization,  this  Hamiltonian  is  identical  in  form  to  the  Hamiltonian  (2)  given  in  the  main 
text.  The  relevant  mechanical  normal  mode  is  given  by  b  =  ( b\  —  b$)/yf2.  Note  that  the  design  is  quite  flexible 
in  that  it  also  applies  to  setups  with  unbalanced  hopping  rates  Kx  ^  K2  and  that,  in  principle,  a  mechanical  mode 
coupled  only  to  either  the  first  or  the  third  site  would  suffice.  The  d0  mode  maybe  driven  through  an  additional 
channel,  provided  that  the  decay  rate  into  this  is  sufficiently  small  so  as  to  ensure  that  the  entangled  photon  pairs 
dominantly  decay  into  another,  outcoupling  waveguide. 

As  pointed  out  in  the  main  text,  a  similar  configuration  of  levels  maybe  realized  in  optomechanical 
membrane-in-the-middle  setups.  In  such  setups,  the  frequencies  of  the  transverse  normal  optical  modes  of  the 
cavity  depend  on  the  membrane  position  and  tilt  angle  [38].  These  parameters  maybe  tuned  so  as  to  create 
triplet  equidistant  optical  modes.  The  frequency  separation  can  be  comparable  to  the  frequency  of  a  vibrational 
mode  of  the  membrane,  which  could  be  matched  to  the  optical  splitting  by  applying  the  optical  spring  effect  [  14] . 
In  such  a  configuration,  the  (linear)  optomechanical  coupling  strengths  are  set  by  the  slopes  of  the  optical  bands 
[37] .  As  pointed  out  above,  even  if  these  turn  out  to  be  different,  that  will  not  impact  our  scheme  in  any 
important  way. 


Appendix  D.  Effective  optical  model 

Next,  we  derive  the  effective  optical  model,  i.e.  the  model  obtained  after  integrating  out  the  mechanics.  We  will 
discuss  how  it  captures  some  essential  features  of  the  full  model.  Assuming  a  large  frequency  mismatch 
6  ft,  T,  g,  |  A  |,  and  low  temperatures,  i.e.,  +  1  <  g/ft,  we  can  adiabatically  eliminate  the  mechanical 

mode  and  the  mechanical  bath.  This  can  be  accomplished  by  a  polaron  transformation  Hopt  =  esHiin  e-S  with 
s  =  2(g  +  A) («+£  -  <U’t)  +  2(6- A) (“-b  -  d-^)- Thus,  we  find 
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Hopt  — 


— A  +  u^ci—)  +  8b  b  —  ^ 


4(5 +A) 

at 


|  (1_|_ 


d2d_  — 


8(5+A) 

.2 


+ 


8(5 -A) 


^(alal  +  d+d_) 


o 


§-a); 


bb 


4(5  -A)"-"  \^4(5+A)  4^' 

;  — A(d+d+  +  dld_)  +  6b' b  —  ^r(d+d+  +  dld_  +  d+d2  +  a+a_). 


(28) 


where,  in  the  last  step,  we  explicitly  used  that  |  A  |  <C  8.  The  optical  modes  are  now  decoupled  from  the 
mechanical  mode  and  we  obtain  a  closed  set  of  quantum  Langevin  equations  for  the  optical  modes 

a+  =  i^A  +  |j)a+  +  i^fll  -  | a+  -  jKa+tin(t), 

«+  =  -i(A  +  S)«+  -  -  f«+  -  VJcat.inW. 

fl_  =  i^A  +  |r)a_  +  -  |fl_  -  VKfl_jin(f), 

fl-  =  -i(a  +  ^)fll  -  i^fl+  -  yfll  -  VKfll  in  (0.  (29) 

where  (aqtin)  =  0,  (aj>in  (t)«9',in  (O)  =  0,  (^.mW^in^')}  =  V^(f  “  Owithq,  q'  =  ±.  The  output 
fields  are  related  by  d±jOUt  (t)  =  -J~k  a±  ( t )  +  d±jin  (t ).  For  notational  convenience,  we  introduce  the  operator 
vectors  A  =  (d+  d_  dl)T  and  Ain  =  (d+jin  d|>in  d_  in  dl  in)T.  The  quantum  Langevin  equations  can  be  cast 
in  the  following  compact  form 

A  =  MA  -  V^Ain,  (30) 

where 


The  system  is  stable  only  if  all  eigenvalues  of  M  have  non-positive  real  parts.  The  boundary  of  stability  is  located 

where  one  of  the  real  parts  becomes  zero,  i.e.,  —  y  +  +  A^A  =  0,  or  equivalently  when 

(A  +  g2/(28))  A  +  k2/ 4  =  0.  This  corresponds  to  the  dashed  line  in  figure  4 (a).  Note  that  the  lower  quadrants 
of  the  stability  diagram  are  essentially  identical  due  to  inversion  symmetry  with  respect  to  the  point  8  =  A  =  0. 

We  can  solve  the  quantum  Langevin  equations  in  Fourier  space.  We  define  the  Fourier-transformed 
operators  by  a+  (cj)  =  f  °°  d+(t)e1LJtdt,  al(— uS)  =  f  °°  dl  (t)elcJtdt  =  In  the  frequency  domain, 

the  input  noise  correlation  reads  (d^in  (cj)d^in(— Cc/))  =  2i r8qq'8(uj  +  u/).  For  convenience,  we  define  the 
vectors  Ain/out  (u)  =  (d+,in/out  (^)d^in/out  (-^)d_,in/out  (w)al>in/out  (-^))T.  The  solution  of  the  Langevin 
equation  can  be  written  as  Aout  (u)  =  S  (u)  Ain  (a;)  with  the  scattering  matrix 

M22  T  'luJ  0  0  — M4  ^ 

0  Mu  +  Mi  4  0 

0  — -Mu  M22  T  i^  0 

Mi 4  0  0  Mu  +  Da 


%)  = 


MV 


(Mi  1  “l-  icj)  (M22  “i-  iw) 


(31) 


where  I  is  the  4x4  identity  matrix.  The  input  noise  correlations  Bin  (u>,  u/)  =  (Ain  (uj)  <g>  A?  (a/) )  are 
scattered  into  Bout  (uj,  uor)  =  (Aout(a;)  (g)  AoTut(^0)  =  S(cj)Bin(u;,  o/)ST(cc/). 

In  the  effective  optical  model,  photons  are  only  created  in  pairs.  The  pair  creation  rate  rpairs  is  equal  to  the 
intensity  (photons  sec-1)  in  any  of  the  two  output  streams.  Due  to  the  choice  of  normalization  in  the  input- 
output  formalism,  this  is  given  by: 


/  1  \  2  p+00  n+oo 

(^±,out^±>out)  =  (yr)  /  du;  /  do;/(d_j_out(  —  cu)d±fOUt 

v  /  —  OO  —  OO 

 (g2/ 46)2k/2 


(k/2)2 


(A  +  £)i 


(32) 
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This  result  shows  that,  depending  on  the  sign  of  ^  A  +  ^  j  A,  we  get  an  enhanced  or  decreased  photon  pair 

creation  rate  compared  with  A  =  0.  In  addition,  since  the  denominator  vanishes  at  the  boundary  of  stability,  the 
photon  pair  creation  rate  diverges  there. 

Appendix  E.  Calculation  of  E  [u] 

Here,  we  review  the  definition  of  the  logarithmic  negativity  and  apply  it  to  quantify  the  entanglement  of  the 
filtered  optical  output  fields,  first  for  the  effective  optical  model  introduced  above,  and  then  for  the  full  model. 

By  applying  narrow  frequency  filters,  we  select  only  single-frequency  components  of  each  of  the  optical  output 
fields.  By  energy  conservation,  correlations  only  occur  between  the  ^-component  of  one  of  the  fields  (say  d+jOUt) 
and  the  —^-component  of  the  other  (d_jOUt).  The  field  operators  for  these  two  single- frequency  output  fields  are 
obtained  as  a+UJ(t)  =  f  f(t  —  s)d+  out  (s)ds,  a _ ^  (t)  =  f*  g(t  —  s)d_  out  (s)ds.  For  convenience,  we 

’  J—oo  ’  ’  oo  °  y _ 

nowwill  choose  Lorentzian  filter  functions  f  (t)  =  ^6  (t)e~^~lujyy  g(t)  =  ^  0  (t)e~^+luJE  with  0(t)  the 
Heaviside  step  function,  as  opposed  to  the  Wannier  basis  used  in  the  main  text.  In  the  limit  of  small  bandwidth 
1  /r  — >  0,  which  is  the  only  one  we  discuss,  however,  both  the  basis  function  probability  [densities  reduce  to 
Dirac  ^-functions,  and  the  results  will  co-incide. 

We  can  use  the  logarithmic  negativity  [9]  to  characterize  the  entanglement  for  the  output  light  beams  [2 1].  In 
order  to  evaluate  it,  we  define  the  vector  u  =  ( x+ ^  ( t)p+  u  (t)x_y_u  ( t)p_  (t))T  with 

Xj(t)  =  +  dj(t)),pj(t)  =  +  (dj(t)  -  aj  (t))(j  =  +,  tv  or  j=  —u>).  The  entanglement  is 

determined  by  the  covariance  matrix  y  with  matrix  elements  =  y( UiUj  +  where  the  operators 

involved  in  this  product  are  all  taken  at  equal  times.  Inserting  the  stationary  solution,  we  find 

‘  n+  0  Re(0  Im(0A 

0  n+  lm(0  — Re(0 

Re(0  lm(0  n_  0 

[lm(0  — Re(0  0  n_ 

where 


V  = 


(33) 


n+  = 


x 

= X 


y—icot  / 


—  oo 
+  oo 


(^+,out(0^+,out(0))  dt  +  1/2  —  |5i4(a;)p  +  1/2, 


(^-,out  (0^-,out  (0))  dt  +  1/2  —  |Si4(o;)|2  +  1/2, 


and 


r 

tJ  —  c 


*  (^+,out  (t)fl— j0ut  (0)  )  dt —  Sn(cj)Si4(  Lj'). 


n+  +  n 


-  )(++)++ being 


The  logarithmic  negativity  is  defined  as  E  =  max  [0,  —In  (2 rj_)  ]  with  rj_  =  2  2 

the  smaller  symplectic  eigenvalue  of  the  partial  transpose  of  matrix  V.  Choosing  u  =  0,  we  plot  E  [u;  =  0]  as  a 
function  of  A/ k  in  figure  4(b). 

To  obtain  the  spectral  density  of  entanglement  in  the  full  model,  we  solve  the  following  system  of  quantum 
Langevin  equations,  which  derive  directly  from  the  model  Hamiltonian  (main  text),  and  where  dissipation  and 
fluctuations  have  been  taken  into  account  using  the  usual  input-output  formalism: 


d+  =  iAd+  +  i  —  yd+  —  Vftd+>in(t), 
fl+  =  -iA al  -  i+  -  y&f.  -  W» 

d_  =  iAd_  +  \^P  —  |d_  —  Vftd_5in(t), 
al  =  —  iAal  —  i  —  yd/  —  VftdjC 
b  =  -iSb  +  i  |(a+  +  al)  —  yb  —  JTbin(t), 


,(0» 


1?  =  i  Sbl  —  i-(«_ 


)  -  T-P  -  VrbX(f), 


(34) 


where,  in  addition,  we  have  {b-m)  =  0,  +  (t)bin  (t1))  =  nthS(t  -  /■')> +  (++'))  =  (nth  +  1)<5  (t  -  t% 
with  =  (exp  (HQ /(k^T))  —  1)_  1 .  For  the  logarithmic  negativity,  we  need  to  evaluate  the  same  optical 
correlations  as  above.  The  results  are  analytical,  but  too  complicated  to  be  reported  here.  Simpler  analytic  results 
can  be  found  for  8  =  A  =  0.  In  this  case  we  find 
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rj_  (ujf  8,  A  =  0)  =  4 C(C  +  +  y)  +  ’  2C^1  +  4^C  +  nth  +  y)  with  the  cooperativity  C  =  yy.The 

entanglement  E  [cj,  <5,  A  =  0]  is  only  determined  by  C  and  nt h  as  given  in  the  main  text. 


Appendix  F.  Comparison  to  two-  and  three-mode  schemes 


Finally,  we  compare  the  efficiency  of  our  four-mode  setup  (3  optical,  1  vibrational)  to  two-  and  three-mode 
schemes  that  have  been  discussed  previously.  We  show  that  the  four-mode  setup  is  more  efficient  than  either  of 
these  schemes. 

For  a  two-mode  setup  (1  optical,  1  vibrational)  and  with  a  resonant  laser  drive  (detuning  A  =  0),  the 
Hamiltonian  is  given  by  H  =  ElPb  —  g(a  +  a^)(b  +  P)>  where  a  is  an  optical  and  bis  a  mechanical  mode 
(see  [14, 21]  for  the  case  of  drive  at  the  red  mechanical  sideband).  The  filtered  correlations  that  quantify  the 
entanglement  between  the  Stokes  and  anti- Stokes  mechanical  sidebands  are  then  given  by 

,  2nth^2 


n  = _ _ , _ i 

(f22  +  (k/2)2)2  (r/2)2  (f22  +  (k/2)2)(T/2)  2’ 


n_  : 


2  4 

Kg 


(O2 


_ 2(»th  +  1  )«g2  i_ 

+  (k/2)2)2  (r/2)2  (O2  +  (k/2)2)  (r/2)  2  ’ 


K2g4 


(2«th  +  1)«£2 


(O2  +  (K/2)2)2(r/2)2  (O2  +  («/2)2)(r/2) 

where  we  have  assumed  that  Q  T.  For  our  model,  with  resonant  drive  A  =  0,  we  find 

nTn^ig/2)1  K2(g/ 2)4 


(35) 


n+  ■■ 


H-  ■ 


£  =  - 


[(-<*>  +  6)2  +  (r/2)2](w2  +  (k/2)2) 

_ «r(»tfa  +  i)(g/2)2 _ 

[(-W  +  6 )2  +  (r/2)2](w2  +  (k/2)2) 

K2(g/2)4 


+  b 


[(-W  +  6)2  +  (r/2)2](w2  +  (k/2)2)2 

_ «2(g/2)4 _ 

[(-LO  +  8)2  +  (r/2)2](w2  +  (k/2)2)2 
+  -Kr(«th  +  l/2)(g/2)2  +  i(— 02  +  8)n(g/2)2 


+  b 


K-U  +  8)2  +  (r/2)2](w2  +  (k/2)2)2  [(-w  +  8 )2  +  (r/2)2](w2  +  (k/2)2) 

In  order  to  evaluate  the  entanglement  between  the  Stokes  and  anti-Stokes  sidebands,  we  set  a;  =  8.  This  yields 

2K«th(g/2)2  K2(g/ 2)4 


(36) 


«+  = 


tl-  = 


£  =  - 


(r/2)(<52  +  (k/2)2) 

2k  (nth  +  l)(g/2)2 

(r/2)(^2  +  (k/2)2) 

K2(g/2)4 


-I - 1 _ l  1 

(r/2)2(<52  +  (k/2)2)2  2’ 

«2(g/2)4  i_ 

(r/2)2(<52  +  (k/2)2)2  2’ 

K(2«th  +  l)(g/2)2 


(r/2)2(<52  +  (k/2)2)2  (r/2)(<52  +  (k/2)2) 

By  comparing  the  coherent  parts  of  n+  and  n_  (the  terms  that  do  not  depend  on  nt h  and  are  useful  for  optical 

Q4 

entanglement),  we  find  that,  in  our  setup,  the  photon  pair  creation  rate  is  enhanced  by  a  factor  ~  +  . 

The  main  benefit  of  our  setup  in  comparison  to  three- mode  schemes,  which  involve  one  mechanical  and 
two  optical  modes,  is  that  the  laser  drive  is  resonant  so  that,  for  a  fixed  drive  strength  A,  the  effective 
optomechanical  coupling  strength  is  much  larger.  By  contrast,  in  a  three- mode  setup,  the  laser  drives  are  off- 
resonant  by  =LQ  so  that  |  a  |  =  ,  A  =,  whereas  in  our  setup,  with  A  =  0,  we  have  |  a  |  =  —.Since  the 


(37) 


b O2  -\-  (k  /  2)" 


intensity  of  entangled  photons  scales  with  g4  and  g  oc  ag0,  this  leads  to  an  enhancement  of  ~  °f  the 

photon  pair  creation  rate. 
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